Overview

Setup

(hide)

Click on tabs to display additional information.

Libraries

# Load required libraries for statistical analysis and table creation
library(dplyr)
library(ggplot2)
library(gt)
library(broom)
library(car)
library(emmeans)
library(multcomp)
library(MASS)  # For negative binomial regression
library(rstatix)
library(coin)
library(phyloseq)
library(microViz)
library(tidyverse)

Plotting

# Define treatment order and color palette
treatment_order <- c(
  "A- T- P-",  # Control
  "A- T- P+",  # Parasite
  "A+ T- P-",  # Antibiotics
  "A+ T- P+",  # Antibiotics_Parasite
  "A- T+ P-",  # Temperature
  "A- T+ P+",  # Temperature_Parasite
  "A+ T+ P-",  # Antibiotics_Temperature
  "A+ T+ P+"   # Antibiotics_Temperature_Parasite
)

# Custom color palette matching treatment order
treatment_colors <- c(
  "#1B9E77",  # A- T- P- (Control)
  "#D95F02",  # A- T- P+ (Parasite)
  "#7570B3",  # A+ T- P- (Antibiotics)
  "#E7298A",  # A+ T- P+ (Antibiotics_Parasite)
  "#66A61E",  # A- T+ P- (Temperature)
  "#E6AB02",  # A- T+ P+ (Temperature_Parasite)
  "#A6761D",  # A+ T+ P- (Antibiotics_Temperature)
  "#666666"   # A+ T+ P+ (Antibiotics_Temperature_Parasite)
)

# Create named vector for color scale
treatment_color_scale <- setNames(treatment_colors, treatment_order)

Functions

# Function to extract sample data as dataframe from phyloseq object
samdatAsDataframe <- function(ps) {
  samdat <- phyloseq::sample_data(ps)
  df <- data.frame(samdat, check.names = FALSE, stringsAsFactors = FALSE)
  return(df)
}

# Function to rename variables in phyloseq object
ps_rename <- function(ps, ...) {
  ps <- microViz::ps_get(ps)
  df <- samdatAsDataframe(ps)
  df <- dplyr::rename(.data = df, ...)
  phyloseq::sample_data(ps) <- df
  return(ps)
}

Import Data

ps.unclean <- readRDS('/Users/michaelsieler/Dropbox/Mac (2)/Documents/Sharpton_Lab/Projects_Repository/Rules_of_Life/major-experiment-2023/Data/Robjects/pseq_uncleaned_05052025.rds')


# Create a dataframe for mortality using the "unclean" phyloseq object
# The unclean phyloseq object is prior to removing samples during post-DADA2 processing (e.g., insufficient reads/sample)

data.mortality <-
    ps.unclean %>%

    ## Update Metadata
    ps_rename(Time = Timepoint) %>%
    microViz::ps_mutate(
        Treatment = case_when(
            Antibiotics == 0 & Temperature == 0 & Pathogen == 0 ~ "A- T- P-",
            Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "A- T- P+",
            Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "A+ T- P-",
            Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "A+ T- P+",
            Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "A- T+ P-",
            Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "A- T+ P+",
            Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "A+ T+ P-",
            Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "A+ T+ P+",
            TRUE ~ "Unknown"
        ), .after = "Pathogen"
    ) %>%
    microViz::ps_mutate(Sample = fecal.sample.number, .before = 1) %>%
    microViz::ps_mutate(Sample = gsub("^f", "", Sample)) %>%
    microViz::ps_filter(Treatment != "Unknown") %>%
    microViz::ps_mutate(
        History = case_when(
            Antibiotics + Temperature == 0 ~ 0,
            Antibiotics + Temperature == 1 ~ 1,
            Antibiotics + Temperature == 2 ~ 2,
        ), .after = "Treatment"
    ) %>%
    
    ## Additional metadata updates, factorizing metadata
    microViz::ps_mutate(
    # Create treatment code
        treatment_code = case_when(
          Antibiotics == 0 & Temperature == 0 & Pathogen == 0 ~ "Aneg_Tneg_Pneg",
          Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "Aneg_Tneg_Ppos",
          Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "Apos_Tneg_Pneg",
          Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "Apos_Tneg_Ppos",
          Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "Aneg_Tpos_Pneg",
          Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "Aneg_Tpos_Ppos",
          Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "Apos_Tpos_Pneg",
          Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "Apos_Tpos_Ppos"
        ),
        # Create treatment group factor
        treatment_group = case_when(
          Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "Parasite",
          Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "Antibiotics",
          Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "Antibiotics_Parasite",
          Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "Temperature",
          Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "Temperature_Parasite",
          Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "Antibiotics_Temperature",
          Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "Antibiotics_Temperature_Parasite",
          TRUE ~ "Control"
        ),
        # Convert to factor with appropriate levels
        treatment_group = factor(treatment_group, 
                               levels = c("Control", "Parasite", 
                                          "Antibiotics", "Antibiotics_Parasite",
                                          "Temperature", "Temperature_Parasite",
                                        "Antibiotics_Temperature", "Antibiotics_Temperature_Parasite")
                               ),
        treatment_code = factor(treatment_code, levels = treatment_order),
        # Create time point factor
        time_point = factor(Time, levels = c(0, 14, 18, 25, 29, 60)),
        # Create pathogen status factor
        pathogen_status = factor(ifelse(Pathogen == 1, "Exposed", "Unexposed"),
                               levels = c("Unexposed", "Exposed")),
        # Create sex factor
        sex = factor(Sex, levels = c("M", "F"))
        ) %>%
    microViz::ps_mutate(Treatment = factor(Treatment, levels = treatment_order)) %>%
      microViz::ps_mutate(Exp_Type = case_when(
          Treatment %in% c("A- T- P-", "A- T- P+")  ~ "No prior stressor(s)",
          Treatment %in% c("A+ T- P-", "A+ T- P+")  ~ "Antibiotics",
          Treatment %in% c("A- T+ P-", "A- T+ P+") ~ "Temperature",
          Treatment %in% c("A+ T+ P-", "A+ T+ P+") ~ "Combined",
      )) %>%
      microViz::ps_mutate(Exp_Type = factor(Exp_Type, levels = c("No prior stressor(s)", "Antibiotics", "Temperature", "Combined"))) %>%
    microViz::samdat_tbl() 


data.infection <- data.mortality %>%
    dplyr::filter(Pathogen == 1) %>%
    dplyr::group_by(Treatment) %>%
    dplyr::summarise(
        n_fish = n(),
        n_infected = sum(Total.Worm.Count > 0),
        mean_worm_burden = mean(Total.Worm.Count, na.rm = TRUE),
        percent_infected = round((n_infected / n_fish) * 100, 1),
        .groups = "drop"
    ) %>%
    dplyr::mutate(
        # Convert Treatment to factor with specific order
        Treatment = factor(Treatment, levels = treatment_order)
    )

data.infection_tank <- data.mortality %>%
    dplyr::filter(Pathogen == 1) %>%
    dplyr::group_by(Treatment, Tank.ID) %>%
    dplyr::summarise(
        n_fish = n(),
        n_infected = sum(Total.Worm.Count > 0),
        mean_worm_burden = mean(Total.Worm.Count, na.rm = TRUE),
        percent_infected = round((n_infected / n_fish) * 100, 1),
        .groups = "drop"
    ) %>%
    dplyr::mutate(
        # Convert Treatment to factor with specific order
        Treatment = factor(Treatment, levels = treatment_order)
    )

# data.mortality
# data.infection
# data.infection_tank

00) All

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality

All

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

By Tank

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

Statistical Analysis - All Treatments Mortality

Code
# Prepare data for statistical analysis - percent mortality rates
mortality_stats_data <- data.mortality %>%
  dplyr::group_by(Treatment) %>%
  dplyr::count() %>%
  dplyr::mutate(
    total_per_group = 90,
    dead = total_per_group - n,
    percent_mortality = round((dead / total_per_group) * 100, 1)
  )

# Chi-square test for independence between treatment and mortality
set.seed(42)
chi_square_test <- chisq.test(
  matrix(c(mortality_stats_data$dead, mortality_stats_data$n), 
         ncol = 2, byrow = FALSE)
)

# Create mortality rate summary table for display
mortality_rate_table <- mortality_stats_data %>%
  dplyr::select(Treatment, dead, n, percent_mortality) %>%
  dplyr::rename(Survived = n) %>%
  dplyr::mutate(
    Total = dead + Survived,
    percent_survived = round((Survived / Total) * 100, 1)
  ) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Mortality Rate Summary - All Treatments",
    subtitle = "Percent mortality and survival rates by treatment"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(dead, Survived, Total),
    decimals = 0
  ) %>%
  gt::fmt_number(
    columns = c(percent_mortality, percent_survived),
    decimals = 1
  )

# Chi-square test results
test_stat_table <- data.frame(
  Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
  Value = c(
    round(chi_square_test$statistic, 4),
    chi_square_test$parameter,
    round(chi_square_test$p.value, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Chi-Square Test Results - All Treatments Mortality",
    subtitle = "Test for independence between treatment and mortality"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Post-hoc pairwise comparisons using chi-square tests
set.seed(42)
pairwise_results <- list()
treatments <- mortality_stats_data$Treatment

for(i in 1:(length(treatments)-1)) {
  for(j in (i+1):length(treatments)) {
    trt1 <- treatments[i]
    trt2 <- treatments[j]
    
    # Get data for these two treatments
    data1 <- mortality_stats_data %>% dplyr::filter(Treatment == trt1)
    data2 <- mortality_stats_data %>% dplyr::filter(Treatment == trt2)
    
    # Create 2x2 contingency table
    cont_table <- matrix(
      c(data1$dead, data1$n, data2$dead, data2$n),
      nrow = 2, ncol = 2,
      dimnames = list(
        c("Dead", "Survived"),
        c(trt1, trt2)
      )
    )
    
    # Perform chi-square test
    test_result <- chisq.test(cont_table)
    
    pairwise_results[[paste(trt1, "vs", trt2)]] <- data.frame(
      Comparison = paste(trt1, "vs", trt2),
      Chi_square = round(test_result$statistic, 4),
      P_value = round(test_result$p.value, 6),
      Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
    )
  }
}

# Combine pairwise results
pairwise_summary <- do.call(rbind, pairwise_results) %>%
    dplyr::arrange(desc(Significant), P_value) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pairwise Chi-Square Tests - All Treatments Mortality",
    subtitle = "Post-hoc comparisons between treatment pairs"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )

# Show only significant pairwise comparisons
pairwise_summary__SigOnly <- do.call(rbind, pairwise_results) %>%
    dplyr::arrange(desc(Significant), P_value) %>%
    dplyr::filter(Significant == "Yes") %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Significant Pairwise Chi-Square Tests - All Treatments Mortality",
    subtitle = "Post-hoc comparisons between treatment pairs (significant only)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )
Table

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality Rates
Mortality Rate Summary - All Treatments
Percent mortality and survival rates by treatment
dead Survived percent_mortality Total percent_survived
A- T- P-
9 81 10.0 90 90.0
A- T- P+
8 82 8.9 90 91.1
A+ T- P-
14 76 15.6 90 84.4
A+ T- P+
10 80 11.1 90 88.9
A- T+ P-
19 71 21.1 90 78.9
A- T+ P+
13 77 14.4 90 85.6
A+ T+ P-
30 60 33.3 90 66.7
A+ T+ P+
20 70 22.2 90 77.8
Chi-Square
Chi-Square Test Results - All Treatments Mortality
Test for independence between treatment and mortality
Statistic Value
Chi-square statistic 29.797800
Degrees of freedom 7.000000
P-value 0.000103
Chi-Square Pairwise
Pairwise Chi-Square Tests - All Treatments Mortality
Post-hoc comparisons between treatment pairs
Comparison Chi_square P_value Significant
A- T- P+ vs A+ T+ P- 14.7109 0.000125 Yes
A- T- P- vs A+ T+ P- 13.0933 0.000296 Yes
A+ T- P+ vs A+ T+ P- 11.6036 0.000658 Yes
A- T+ P+ vs A+ T+ P- 7.8221 0.005161 Yes
A+ T- P- vs A+ T+ P- 6.7680 0.009280 Yes
A- T- P+ vs A+ T+ P+ 5.1175 0.023686 Yes
A- T- P+ vs A- T+ P- 4.3573 0.036851 Yes
A- T- P- vs A+ T+ P+ 4.1105 0.042617 Yes
A- T- P- vs A- T+ P- 3.4258 0.064187 No
A+ T- P+ vs A+ T+ P+ 3.2400 0.071861 No
A- T+ P- vs A+ T+ P- 2.8042 0.094019 No
A+ T- P+ vs A- T+ P- 2.6307 0.104813 No
A+ T+ P- vs A+ T+ P+ 2.2431 0.134214 No
A- T+ P+ vs A+ T+ P+ 1.3358 0.247775 No
A- T- P+ vs A+ T- P- 1.2946 0.255204 No
A- T+ P- vs A- T+ P+ 0.9502 0.329676 No
A+ T- P- vs A+ T+ P+ 0.9065 0.341038 No
A- T- P+ vs A- T+ P+ 0.8625 0.353031 No
A- T- P- vs A+ T- P- 0.7976 0.371823 No
A+ T- P- vs A- T+ P- 0.5937 0.440995 No
A- T- P- vs A- T+ P+ 0.4661 0.494809 No
A+ T- P- vs A+ T- P+ 0.4327 0.510671 No
A+ T- P+ vs A- T+ P+ 0.1994 0.655213 No
A- T- P+ vs A+ T- P+ 0.0617 0.803785 No
A- T- P- vs A- T- P+ 0.0000 1.000000 No
A- T- P- vs A+ T- P+ 0.0000 1.000000 No
A+ T- P- vs A- T+ P+ 0.0000 1.000000 No
A- T+ P- vs A+ T+ P+ 0.0000 1.000000 No
Significant Pairwise Chi-Square Tests - All Treatments Mortality
Post-hoc comparisons between treatment pairs (significant only)
Comparison Chi_square P_value Significant
A- T- P+ vs A+ T+ P- 14.7109 0.000125 Yes
A- T- P- vs A+ T+ P- 13.0933 0.000296 Yes
A+ T- P+ vs A+ T+ P- 11.6036 0.000658 Yes
A- T+ P+ vs A+ T+ P- 7.8221 0.005161 Yes
A+ T- P- vs A+ T+ P- 6.7680 0.009280 Yes
A- T- P+ vs A+ T+ P+ 5.1175 0.023686 Yes
A- T- P+ vs A- T+ P- 4.3573 0.036851 Yes
A- T- P- vs A+ T+ P+ 4.1105 0.042617 Yes

01) What is the host-microbiome response to parasite exposure?

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality

All

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

By Tank

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

Statistical Analysis - Parasite Exposure Mortality

Code
# Prepare data for statistical analysis - parasite exposure comparison (percent mortality rates)
parasite_mortality_data <- data.mortality %>%
  dplyr::filter(Treatment %in% c("A- T- P-", "A- T- P+")) %>%
  dplyr::group_by(Treatment) %>%
  dplyr::count() %>%
  dplyr::mutate(
    total_per_group = 90,
    dead = total_per_group - n,
    percent_mortality = round((dead / total_per_group) * 100, 1)
  )

# Fisher's exact test for 2x2 contingency table (more appropriate for small samples)
set.seed(42)
fisher_test <- fisher.test(
  matrix(c(parasite_mortality_data$dead, parasite_mortality_data$n), 
         ncol = 2, byrow = FALSE)
)

# Chi-square test as alternative
chi_square_parasite <- chisq.test(
  matrix(c(parasite_mortality_data$dead, parasite_mortality_data$n), 
         ncol = 2, byrow = FALSE)
)

# Create mortality rate summary table for display
parasite_mortality_rate_table <- parasite_mortality_data %>%
  dplyr::select(Treatment, dead, n, percent_mortality) %>%
  dplyr::rename(Survived = n) %>%
  dplyr::mutate(
    Total = dead + Survived,
    percent_survived = round((Survived / Total) * 100, 1)
  ) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Mortality Rate Summary - Parasite Exposure",
    subtitle = "Percent mortality and survival rates: A-T-P- vs A-T-P+"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(dead, Survived, Total),
    decimals = 0
  ) %>%
  gt::fmt_number(
    columns = c(percent_mortality, percent_survived),
    decimals = 1
  )

# Fisher's exact test results
fisher_results <- data.frame(
  Statistic = c("Fisher's exact test p-value", "Odds ratio", "95% CI lower", "95% CI upper"),
  Value = c(
    round(fisher_test$p.value, 6),
    round(fisher_test$estimate, 4),
    round(fisher_test$conf.int[1], 4),
    round(fisher_test$conf.int[2], 4)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Fisher's Exact Test Results - Parasite Exposure Mortality",
    subtitle = "Comparison of A-T-P- vs A-T-P+ mortality rates"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Chi-square test results for comparison
chi_square_parasite_results <- data.frame(
  Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
  Value = c(
    round(chi_square_parasite$statistic, 4),
    chi_square_parasite$parameter,
    round(chi_square_parasite$p.value, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Chi-Square Test Results - Parasite Exposure Mortality",
    subtitle = "Alternative test for A-T-P- vs A-T-P+ mortality rates"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Additional analysis: Compare mortality rates between parasite-exposed and unexposed
# Calculate difference in mortality rates
parasite_mortality_diff <- parasite_mortality_data %>%
  dplyr::arrange(Treatment) %>%
  dplyr::mutate(
    Pathogen_Status = dplyr::case_when(
      stringr::str_detect(Treatment, "P\\+") ~ "Exposed",
      stringr::str_detect(Treatment, "P\\-") ~ "Unexposed",
      TRUE ~ "Unknown"
    )
  ) %>%
  dplyr::select(Pathogen_Status, percent_mortality)
## Adding missing grouping variables: `Treatment`
# Calculate the difference in mortality rates
mortality_rate_diff <- parasite_mortality_diff$percent_mortality[2] - parasite_mortality_diff$percent_mortality[1]

# Create summary of mortality rate difference
parasite_mortality_diff_summary <- data.frame(
  Comparison = "A-T-P+ vs A-T-P-",
  Mortality_Rate_Difference = round(mortality_rate_diff, 2),
  Unexposed_Rate = round(parasite_mortality_diff$percent_mortality[1], 2),
  Exposed_Rate = round(parasite_mortality_diff$percent_mortality[2], 2),
  Percent_Increase = round((mortality_rate_diff / parasite_mortality_diff$percent_mortality[1]) * 100, 1)
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Mortality Rate Difference - Parasite Exposure",
    subtitle = "Comparison of mortality rates between parasite-exposed and unexposed groups"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(Mortality_Rate_Difference, Unexposed_Rate, Exposed_Rate, Percent_Increase),
    decimals = 2
  )
Table

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality Rates
Mortality Rate Summary - Parasite Exposure
Percent mortality and survival rates: A-T-P- vs A-T-P+
dead Survived percent_mortality Total percent_survived
A- T- P-
9 81 10.0 90 90.0
A- T- P+
8 82 8.9 90 91.1
Mortality Rate Difference
Mortality Rate Difference - Parasite Exposure
Comparison of mortality rates between parasite-exposed and unexposed groups
Comparison Mortality_Rate_Difference Unexposed_Rate Exposed_Rate Percent_Increase
A-T-P+ vs A-T-P- −1.10 10.00 8.90 −11.00
Fisher’s Exact
Fisher's Exact Test Results - Parasite Exposure Mortality
Comparison of A-T-P- vs A-T-P+ mortality rates
Statistic Value
Fisher's exact test p-value 1.0000
Odds ratio 1.1381
95% CI lower 0.3692
95% CI upper 3.5739
Chi-Square
Chi-Square Test Results - Parasite Exposure Mortality
Alternative test for A-T-P- vs A-T-P+ mortality rates
Statistic Value
Chi-square statistic 0
Degrees of freedom 1
P-value 1

02) Is the host-microbiome response to parasite exposure historically contingent?

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality

All

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

By Tank

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

Statistical Analysis - Historical Contingency Mortality

Code
# Prepare data for statistical analysis - historical contingency mortality
historical_mortality_data <- data.mortality %>%
  dplyr::filter(Treatment %in% c("A- T- P+", "A+ T- P+", "A- T+ P+", "A+ T+ P+")) %>%
  dplyr::group_by(Treatment) %>%
  dplyr::count() %>%
  dplyr::mutate(
    total_per_group = 90,
    dead = total_per_group - n,
    percent_mortality = round((dead / total_per_group) * 100, 1)
  )

# Chi-square test for independence between treatment and mortality
set.seed(42)
historical_chi_square_test <- chisq.test(
  matrix(c(historical_mortality_data$dead, historical_mortality_data$n), 
         ncol = 2, byrow = FALSE)
)

# Create mortality rate summary table for display
historical_mortality_rate_table <- historical_mortality_data %>%
  dplyr::select(Treatment, dead, n, percent_mortality) %>%
  dplyr::rename(Survived = n) %>%
  dplyr::mutate(
    Total = dead + Survived,
    percent_survived = round((Survived / Total) * 100, 1)
  ) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Mortality Rate Summary - Historical Contingency",
    subtitle = "Percent mortality and survival rates by treatment (parasite-exposed only)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(dead, Survived, Total),
    decimals = 0
  ) %>%
  gt::fmt_number(
    columns = c(percent_mortality, percent_survived),
    decimals = 1
  )

# Chi-square test results
historical_chi_square_results <- data.frame(
  Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
  Value = c(
    round(historical_chi_square_test$statistic, 4),
    historical_chi_square_test$parameter,
    round(historical_chi_square_test$p.value, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Chi-Square Test Results - Historical Contingency Mortality",
    subtitle = "Test for independence between treatment and mortality (parasite-exposed only)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Post-hoc pairwise comparisons using chi-square tests
set.seed(42)
historical_pairwise_results <- list()
historical_treatments <- historical_mortality_data$Treatment

for(i in 1:(length(historical_treatments)-1)) {
  for(j in (i+1):length(historical_treatments)) {
    trt1 <- historical_treatments[i]
    trt2 <- historical_treatments[j]
    
    # Get data for these two treatments
    data1 <- historical_mortality_data %>% dplyr::filter(Treatment == trt1)
    data2 <- historical_mortality_data %>% dplyr::filter(Treatment == trt2)
    
    # Create 2x2 contingency table
    cont_table <- matrix(
      c(data1$dead, data1$n, data2$dead, data2$n),
      nrow = 2, ncol = 2,
      dimnames = list(
        c("Dead", "Survived"),
        c(trt1, trt2)
      )
    )
    
    # Perform chi-square test
    test_result <- chisq.test(cont_table)
    
    historical_pairwise_results[[paste(trt1, "vs", trt2)]] <- data.frame(
      Comparison = paste(trt1, "vs", trt2),
      Chi_square = round(test_result$statistic, 4),
      P_value = round(test_result$p.value, 6),
      Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
    )
  }
}

# Combine pairwise results
historical_pairwise_summary <- do.call(rbind, historical_pairwise_results) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pairwise Chi-Square Tests - Historical Contingency Mortality",
    subtitle = "Post-hoc comparisons between treatment pairs (parasite-exposed only)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )


data.mortality %>%
    dplyr::filter(Treatment %in% c("A- T- P+", "A+ T- P+", "A- T+ P+", "A+ T+ P+")) %>%
  dplyr::group_by(Treatment) %>%
  dplyr::count() %>%
  dplyr::mutate(
    total_per_group = 90,
    dead = total_per_group - n,
    percent_mortality = round((dead / total_per_group) * 100, 1)
  )
## # A tibble: 4 × 5
## # Groups:   Treatment [4]
##   Treatment     n total_per_group  dead percent_mortality
##   <fct>     <int>           <dbl> <dbl>             <dbl>
## 1 A- T- P+     82              90     8               8.9
## 2 A+ T- P+     80              90    10              11.1
## 3 A- T+ P+     77              90    13              14.4
## 4 A+ T+ P+     70              90    20              22.2
Table

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality Rates
Mortality Rate Summary - Historical Contingency
Percent mortality and survival rates by treatment (parasite-exposed only)
dead Survived percent_mortality Total percent_survived
A- T- P+
8 82 8.9 90 91.1
A+ T- P+
10 80 11.1 90 88.9
A- T+ P+
13 77 14.4 90 85.6
A+ T+ P+
20 70 22.2 90 77.8
Chi-Square
Chi-Square Test Results - Historical Contingency Mortality
Test for independence between treatment and mortality (parasite-exposed only)
Statistic Value
Chi-square statistic 7.561400
Degrees of freedom 3.000000
P-value 0.056002
Chi-Square Pairwise
Pairwise Chi-Square Tests - Historical Contingency Mortality
Post-hoc comparisons between treatment pairs (parasite-exposed only)
Comparison Chi_square P_value Significant
A- T- P+ vs A+ T- P+ 0.0617 0.803785 No
A- T- P+ vs A- T+ P+ 0.8625 0.353031 No
A- T- P+ vs A+ T+ P+ 5.1175 0.023686 Yes
A+ T- P+ vs A- T+ P+ 0.1994 0.655213 No
A+ T- P+ vs A+ T+ P+ 3.2400 0.071861 No
A- T+ P+ vs A+ T+ P+ 1.3358 0.247775 No

Infection

All

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

By Tank

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

Statistical Analysis - Infection Prevalence and Worm Burden

Code
# Statistical analysis for infection prevalence and worm burden
set.seed(42)

# 1. INFECTION PREVALENCE ANALYSIS
# Chi-square test for infection prevalence across treatments
infection_chi_square <- chisq.test(
  data.infection %>%
    dplyr::select(n_infected, n_fish) %>%
    dplyr::mutate(n_uninfected = n_fish - n_infected) %>%
    dplyr::select(n_infected, n_uninfected) %>%
    as.matrix()
)

# Create infection prevalence contingency table
infection_contingency_table <- data.infection %>%
  dplyr::mutate(
    n_uninfected = n_fish - n_infected,
    percent_uninfected = round((n_uninfected / n_fish) * 100, 1)
  ) %>%
  dplyr::select(Treatment, n_infected, n_uninfected, n_fish, percent_infected, percent_uninfected) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Infection Prevalence Contingency Table",
    subtitle = "Number of infected vs uninfected fish by treatment"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(n_infected, n_uninfected, n_fish),
    decimals = 0
  )

# Chi-square test results for infection prevalence
infection_chi_results <- data.frame(
  Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
  Value = c(
    round(infection_chi_square$statistic, 4),
    infection_chi_square$parameter,
    round(infection_chi_square$p.value, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Chi-Square Test Results - Infection Prevalence",
    subtitle = "Test for independence between treatment and infection status"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Post-hoc pairwise chi-square tests for infection prevalence
set.seed(42)
infection_pairwise_results <- list()
infection_treatments <- data.infection$Treatment

for(i in 1:(length(infection_treatments)-1)) {
  for(j in (i+1):length(infection_treatments)) {
    trt1 <- infection_treatments[i]
    trt2 <- infection_treatments[j]
    
    # Get data for these two treatments
    data1 <- data.infection %>% dplyr::filter(Treatment == trt1)
    data2 <- data.infection %>% dplyr::filter(Treatment == trt2)
    
    # Create 2x2 contingency table
    cont_table <- matrix(
      c(data1$n_infected, data1$n_fish - data1$n_infected, 
        data2$n_infected, data2$n_fish - data2$n_infected),
      nrow = 2, ncol = 2,
      dimnames = list(
        c("Infected", "Uninfected"),
        c(trt1, trt2)
      )
    )
    
    # Perform chi-square test
    test_result <- chisq.test(cont_table)
    
    infection_pairwise_results[[paste(trt1, "vs", trt2)]] <- data.frame(
      Comparison = paste(trt1, "vs", trt2),
      Chi_square = round(test_result$statistic, 4),
      P_value = round(test_result$p.value, 6),
      Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
    )
  }
}

# Combine pairwise results for infection prevalence
infection_pairwise_summary <- do.call(rbind, infection_pairwise_results) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pairwise Chi-Square Tests - Infection Prevalence",
    subtitle = "Post-hoc comparisons between treatment pairs"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )

# 2. WORM BURDEN ANALYSIS (TOTAL WORMS)
# Prepare individual-level data for negative binomial regression
set.seed(42)
individual_worm_data <- data.mortality %>%
  dplyr::filter(Pathogen == 1) %>%
  dplyr::select(Treatment, Total.Worm.Count) %>%
  dplyr::filter(!is.na(Total.Worm.Count)) %>%
  dplyr::mutate(
    Treatment = factor(Treatment, levels = c("A- T- P+", "A+ T- P+", "A- T+ P+", "A+ T+ P+"))
  )

# Fit negative binomial model for total worm counts
# Using MASS::glm.nb for negative binomial regression
library(MASS)
nb_model <- MASS::glm.nb(Total.Worm.Count ~ Treatment, data = individual_worm_data)

# Get model summary
nb_summary <- summary(nb_model)

nb_model.anova <- Anova(nb_model, type = 2)

# Extract coefficients and significance
nb_coefficients <- data.frame(
  Treatment = c("Intercept", levels(individual_worm_data$Treatment)[-1]),
  Estimate = round(nb_summary$coefficients[, 1], 4),
  Std_Error = round(nb_summary$coefficients[, 2], 4),
  z_value = round(nb_summary$coefficients[, 3], 4),
  p_value = round(nb_summary$coefficients[, 4], 6)
) %>%
  dplyr::mutate(
    p_signif = ifelse(p_value < 0.001, "***", 
                     ifelse(p_value < 0.01, "**",
                            ifelse(p_value < 0.05, "*", "ns")))
  )

# Create worm burden summary table
worm_burden_data <- data.mortality %>%
  dplyr::filter(Pathogen == 1) %>%
  dplyr::group_by(Treatment) %>%
  dplyr::summarise(
    total_worms = sum(Total.Worm.Count, na.rm = TRUE),
    mean_worms = mean(Total.Worm.Count, na.rm = TRUE),
    median_worms = median(Total.Worm.Count, na.rm = TRUE),
    sd_worms = sd(Total.Worm.Count, na.rm = TRUE),
    n_fish = n(),
    .groups = "drop"
  ) %>%
  dplyr::mutate(
    se_worms = sd_worms / sqrt(n_fish)
  )

worm_burden_summary <- worm_burden_data %>%
  dplyr::select(Treatment, total_worms, mean_worms, median_worms, sd_worms, se_worms, n_fish) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Worm Burden Summary by Treatment",
    subtitle = "Total, mean, median, standard deviation, and standard error of worm counts per treatment"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(total_worms, n_fish),
    decimals = 0
  ) %>%
  gt::fmt_number(
    columns = c(mean_worms, median_worms, sd_worms, se_worms),
    decimals = 2
  )

# Negative binomial model results
nb_model_results <- nb_coefficients %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Negative Binomial Regression Results - Worm Burden",
    subtitle = "Model: glm.nb(worm counts ~ Treatment)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(Estimate, Std_Error, z_value, p_value),
    decimals = 3
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = p_signif,
      rows = p_signif != "ns"
    )
  )


nb_model.anova %>%
    broom::tidy() %>%
  dplyr::mutate(
    p_signif = ifelse(p.value < 0.001, "***", 
                     ifelse(p.value < 0.01, "**",
                            ifelse(p.value < 0.05, "*", "ns")))
  ) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Negative Binomial Regression ANOVA Results - Worm Burden",
    subtitle = "Model: Anova(glm.nb(worm counts ~ Treatment))"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(p.value),
    decimals = 3
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = p_signif,
      rows = p_signif != "ns"
    )
  )
Negative Binomial Regression ANOVA Results - Worm Burden
Model: Anova(glm.nb(worm counts ~ Treatment))
term statistic df p.value p_signif
Treatment 9.044622 3 0.029 *
# Model fit statistics
nb_fit_stats <- data.frame(
  Statistic = c("Log-likelihood", "AIC", "Theta", "SE Theta", "Deviance", "Residual Deviance"),
  Value = c(
    round(logLik(nb_model), 2),
    round(AIC(nb_model), 2),
    round(nb_model$theta, 4),
    round(nb_model$SE.theta, 4),
    round(nb_summary$deviance, 2),
    round(sum(residuals(nb_model, type = "deviance")^2), 2)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Negative Binomial Model Fit Statistics",
    subtitle = "Model diagnostics and fit measures"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Post-hoc pairwise comparisons for worm burden using emmeans
set.seed(42)
library(emmeans)

# Get estimated marginal means and pairwise comparisons
nb_emmeans <- emmeans::emmeans(nb_model, ~ Treatment)
nb_pairwise <- pairs(nb_emmeans, adjust = "bonferroni")

# Convert to data frame for display
nb_pairwise_df <- as.data.frame(nb_pairwise) %>%
  dplyr::mutate(
    contrast = as.character(contrast),
    estimate = round(estimate, 4),
    SE = round(SE, 4),
    z.ratio = round(z.ratio, 4),
    p.value = round(p.value, 6),
    p_signif = ifelse(p.value < 0.001, "***", 
                     ifelse(p.value < 0.01, "**",
                            ifelse(p.value < 0.05, "*", "ns")))
  ) %>%
  dplyr::select(contrast, estimate, SE, z.ratio, p.value, p_signif)

# Display pairwise comparisons for worm burden
nb_pairwise_summary <- nb_pairwise_df %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pairwise Comparisons - Worm Burden (Negative Binomial)",
    subtitle = "Bonferroni-adjusted pairwise comparisons between treatments"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(estimate, SE, z.ratio, p.value),
    decimals = 3
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = p_signif,
      rows = p_signif != "ns"
    )
  )

# Model diagnostics for negative binomial
# Check dispersion parameter and model adequacy
nb_diagnostics <- data.frame(
  Statistic = c("Theta (dispersion parameter)", "SE Theta", "Deviance", "Degrees of Freedom", "Deviance/DF", "Model adequacy"),
  Value = c(
    round(nb_model$theta, 4),
    round(nb_model$SE.theta, 4),
    round(nb_summary$deviance, 2),
    nb_summary$df.residual,
    round(nb_summary$deviance / nb_summary$df.residual, 2),
    ifelse(nb_summary$deviance / nb_summary$df.residual < 2, "Good fit", "Consider model diagnostics")
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Negative Binomial Model Diagnostics",
    subtitle = "Checking model adequacy and dispersion"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Compare with Poisson model for reference
set.seed(42)
poisson_model <- glm(Total.Worm.Count ~ Treatment, data = individual_worm_data, family = "poisson")
poisson_summary <- summary(poisson_model)

# Model comparison
model_comparison <- data.frame(
  Model = c("Poisson", "Negative Binomial"),
  Log_Likelihood = c(round(logLik(poisson_model), 2), round(logLik(nb_model), 2)),
  AIC = c(round(AIC(poisson_model), 2), round(AIC(nb_model), 2)),
  Deviance = c(round(poisson_summary$deviance, 2), round(nb_summary$deviance, 2)),
  DF = c(poisson_summary$df.residual, nb_summary$df.residual),
  Deviance_DF_Ratio = c(
    round(poisson_summary$deviance / poisson_summary$df.residual, 2),
    round(nb_summary$deviance / nb_summary$df.residual, 2)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Model Comparison: Poisson vs Negative Binomial",
    subtitle = "Comparing model fit and adequacy"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(Log_Likelihood, AIC, Deviance, Deviance_DF_Ratio),
    decimals = 2
  ) %>%
  gt::fmt_number(
    columns = DF,
    decimals = 0
  )

# 3. NORMALIZED WORM BURDEN ANALYSIS (WORMS PER INFECTED FISH)
# Calculate worms per infected fish to normalize for infection prevalence
# NOTE: This is ratio data, not count data, so we need different statistical approaches
set.seed(42)

# Create normalized worm burden summary table
normalized_worm_summary <- normalized_worm_data %>%
  dplyr::select(Treatment, total_worms, n_infected, worms_per_infected, n_fish) %>%
  dplyr::mutate(
    infection_rate = round((n_infected / n_fish) * 100, 1)
  ) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Normalized Worm Burden Summary by Treatment",
    subtitle = "Worms per infected fish (normalized for infection prevalence)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(total_worms, n_infected, n_fish),
    decimals = 0
  ) %>%
  gt::fmt_number(
    columns = c(worms_per_infected, infection_rate),
    decimals = 2
  )

# Create tank_normalized_for_test dataset for statistical analysis
# Using ANOVA on worms per infected fish (excluding tanks with 0 infected fish)
set.seed(42)
tank_normalized_for_test <- tank_normalized_data %>%
  dplyr::filter(n_infected > 0) %>%
  dplyr::mutate(
    Treatment = factor(Treatment, levels = treatment_order)
  )

# APPROPRIATE STATISTICAL ANALYSIS FOR RATIO DATA
# Method 1: One-way ANOVA on tank-level ratios (if assumptions are met)
set.seed(42)
# Use tank-level data for statistical testing
normalized_anova <- aov(worms_per_infected ~ Treatment, data = tank_normalized_for_test)
normalized_anova_summary <- summary(normalized_anova)

# Extract ANOVA results using broom package for cleaner extraction
library(broom)
normalized_anova_tidy <- tidy(normalized_anova)

# Calculate F-statistic and p-value manually from ANOVA summary
f_stat <- normalized_anova_summary[[1]]$"Mean Sq"[1] / normalized_anova_summary[[1]]$"Mean Sq"[2]
p_val <- pf(f_stat, normalized_anova_summary[[1]]$"Df"[1], normalized_anova_summary[[1]]$"Df"[2], lower.tail = FALSE)

normalized_anova_results <- data.frame(
  Statistic = c("F-value", "Degrees of freedom (Treatment)", "Degrees of freedom (Residuals)", "P-value"),
  Value = c(
    round(f_stat, 4),
    normalized_anova_summary[[1]]$"Df"[1],
    normalized_anova_summary[[1]]$"Df"[2],
    round(p_val, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "ANOVA Results - Normalized Worm Burden",
    subtitle = "One-way ANOVA: worms per infected fish ~ Treatment (tank-level data)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Post-hoc Tukey test for normalized worm burden
set.seed(42)
normalized_tukey <- TukeyHSD(normalized_anova)

# Convert Tukey results to data frame
normalized_tukey_df <- as.data.frame(normalized_tukey$Treatment) %>%
  dplyr::mutate(
    comparison = rownames(.),
    diff = round(diff, 4),
    lwr = round(lwr, 4),
    upr = round(upr, 4),
    p_adj = round(`p adj`, 6),
    p_signif = ifelse(p_adj < 0.001, "***",
                     ifelse(p_adj < 0.01, "**",
                            ifelse(p_adj < 0.05, "*", "ns")))
  ) %>%
  dplyr::select(comparison, diff, lwr, upr, p_adj, p_signif)

# Display Tukey test results
normalized_tukey_summary <- normalized_tukey_df %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Tukey HSD Test Results - Normalized Worm Burden",
    subtitle = "Post-hoc pairwise comparisons for worms per infected fish (tank-level data)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(diff, lwr, upr, p_adj),
    decimals = 3
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = p_signif,
      rows = p_signif != "ns"
    )
  )

# Method 2: Kruskal-Wallis test (non-parametric alternative)
set.seed(42)
normalized_kruskal <- kruskal.test(worms_per_infected ~ Treatment, data = tank_normalized_for_test)

# Kruskal-Wallis results
normalized_kruskal_results <- data.frame(
  Statistic = c("Kruskal-Wallis chi-squared", "Degrees of freedom", "P-value"),
  Value = c(
    round(normalized_kruskal$statistic, 4),
    normalized_kruskal$parameter,
    round(normalized_kruskal$p.value, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Kruskal-Wallis Test Results - Normalized Worm Burden",
    subtitle = "Non-parametric test for differences in worms per infected fish across treatments (tank-level data)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Method 3: Permutation test (most appropriate for ratio data with small sample sizes)
set.seed(42)
library(coin)

# Create permutation test for normalized worm burden using tank-level data
normalized_permutation <- coin::oneway_test(
  worms_per_infected ~ Treatment,
  data = tank_normalized_for_test,
  distribution = approximate(nresample = 10000)
)

# Permutation test results
normalized_permutation_results <- data.frame(
  Statistic = c("Z statistic", "P-value (approximate)", "Number of permutations"),
  Value = c(
    round(statistic(normalized_permutation), 4),
    round(pvalue(normalized_permutation), 6),
    attr(pvalue(normalized_permutation), "nresample")
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Permutation Test Results - Normalized Worm Burden",
    subtitle = "Non-parametric permutation test for ratio data (tank-level data)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Method 4: Bootstrap confidence intervals for ratios
set.seed(42)
library(boot)
## 
## Attaching package: 'boot'
## The following object is masked from 'package:survival':
## 
##     aml
## The following object is masked from 'package:car':
## 
##     logit
# Function to calculate mean ratio for bootstrap
ratio_mean <- function(data, indices) {
  d <- data[indices, ]
  return(mean(d$worms_per_infected, na.rm = TRUE))
}

# Bootstrap analysis for each treatment using tank-level data
# We need multiple observations per treatment for bootstrap to work
bootstrap_results <- list()
for(trt in unique(tank_normalized_data$Treatment)) {
  trt_data <- tank_normalized_data %>% dplyr::filter(Treatment == trt, n_infected > 0)

  # Only proceed if we have multiple observations
  if(nrow(trt_data) > 1) {
    # Bootstrap with 10000 resamples
    boot_result <- boot(trt_data, ratio_mean, R = 10000)

    # Calculate confidence intervals
    ci_result <- boot.ci(boot_result, type = "perc", conf = 0.95)

    bootstrap_results[[trt]] <- data.frame(
      Treatment = trt,
      Mean_Ratio = round(mean(trt_data$worms_per_infected, na.rm = TRUE), 3),
      Bootstrap_Mean = round(boot_result$t0, 3),
      CI_Lower = round(ci_result$percent[4], 3),
      CI_Upper = round(ci_result$percent[5], 3),
      N_Tanks = nrow(trt_data)
    )
  } else {
    # For treatments with only one observation, just report the mean
    bootstrap_results[[trt]] <- data.frame(
      Treatment = trt,
      Mean_Ratio = round(mean(trt_data$worms_per_infected, na.rm = TRUE), 3),
      Bootstrap_Mean = round(mean(trt_data$worms_per_infected, na.rm = TRUE), 3),
      CI_Lower = NA,
      CI_Upper = NA,
      N_Tanks = nrow(trt_data)
    )
  }
}

# Combine bootstrap results
bootstrap_summary <- do.call(rbind, bootstrap_results) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Bootstrap Confidence Intervals - Normalized Worm Burden",
    subtitle = "95% confidence intervals for worms per infected fish by treatment (tank-level data)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(Mean_Ratio, Bootstrap_Mean, CI_Lower, CI_Upper),
    decimals = 3
  )

# Method 5: Pairwise permutation tests
set.seed(42)
pairwise_permutation_results <- list()
treatments <- unique(tank_normalized_for_test$Treatment)

for(i in 1:(length(treatments)-1)) {
  for(j in (i+1):length(treatments)) {
    trt1 <- treatments[i]
    trt2 <- treatments[j]

    # Get data for these two treatments
    data1 <- tank_normalized_for_test %>% dplyr::filter(Treatment == trt1)
    data2 <- tank_normalized_for_test %>% dplyr::filter(Treatment == trt2)

    # Combine data for permutation test
    combined_data <- rbind(data1, data2)

    # Perform permutation test
    perm_test <- coin::oneway_test(
      worms_per_infected ~ Treatment,
      data = combined_data,
      distribution = approximate(nresample = 10000)
    )

    pairwise_permutation_results[[paste(trt1, "vs", trt2)]] <- data.frame(
      Comparison = paste(trt1, "vs", trt2),
      Z_statistic = round(statistic(perm_test), 4),
      P_value = round(pvalue(perm_test), 6),
      Significant = ifelse(pvalue(perm_test) < 0.05, "Yes", "No")
    )
  }
}

# Combine pairwise permutation results
pairwise_permutation_summary <- do.call(rbind, pairwise_permutation_results) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pairwise Permutation Tests - Normalized Worm Burden",
    subtitle = "Non-parametric pairwise comparisons for ratio data (tank-level data)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )

# Model comparison and diagnostics
normalized_model_comparison <- data.frame(
  Test = c("ANOVA", "Kruskal-Wallis", "Permutation Test"),
  Statistic = c(
    paste("F =", round(f_stat, 4)),
    paste("χ² =", round(normalized_kruskal$statistic, 4)),
    paste("Z =", round(statistic(normalized_permutation), 4))
  ),
  P_value = c(
    round(p_val, 6),
    round(normalized_kruskal$p.value, 6),
    round(pvalue(normalized_permutation), 6)
  ),
  Assumptions = c(
    "Normality, homogeneity of variance",
    "None (non-parametric)",
    "None (non-parametric)"
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Model Comparison - Normalized Worm Burden Analysis",
    subtitle = "Comparing different statistical approaches for ratio data"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Tank-level normalized worm burden analysis
set.seed(42)

# Tank-level normalized worm burden summary
tank_normalized_summary <- tank_normalized_data %>%
  dplyr::select(Treatment, Tank.ID, total_worms, n_infected, worms_per_infected, n_fish, infection_rate) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Tank-Level Normalized Worm Burden Summary",
    subtitle = "Worms per infected fish by treatment and tank"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(total_worms, n_infected, n_fish),
    decimals = 0
  ) %>%
  gt::fmt_number(
    columns = c(worms_per_infected, infection_rate),
    decimals = 2
  )

# Statistical test for tank-level normalized worm burden
# Using ANOVA on worms per infected fish (excluding tanks with 0 infected fish)
set.seed(42)

# One-way ANOVA for worms per infected fish across treatments
tank_normalized_anova <- aov(worms_per_infected ~ Treatment, data = tank_normalized_for_test)
tank_normalized_anova_summary <- summary(tank_normalized_anova)

# Extract ANOVA results
tank_normalized_anova_results <- data.frame(
  Statistic = c("F-value", "Degrees of freedom (Treatment)", "Degrees of freedom (Residuals)", "P-value"),
  Value = c(
    round(tank_normalized_anova_summary[[1]]$"F value"[1], 4),
    tank_normalized_anova_summary[[1]]$"Df"[1],
    tank_normalized_anova_summary[[1]]$"Df"[2],
    round(tank_normalized_anova_summary[[1]]$"Pr(>F)"[1], 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "ANOVA Results - Tank-Level Normalized Worm Burden",
    subtitle = "One-way ANOVA: worms per infected fish ~ Treatment"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Post-hoc Tukey test for tank-level normalized worm burden
set.seed(42)
tank_normalized_tukey <- TukeyHSD(tank_normalized_anova)

# Convert Tukey results to data frame
tank_normalized_tukey_df <- as.data.frame(tank_normalized_tukey$Treatment) %>%
  dplyr::mutate(
    comparison = rownames(.),
    diff = round(diff, 4),
    lwr = round(lwr, 4),
    upr = round(upr, 4),
    p_adj = round(`p adj`, 6),
    p_signif = ifelse(p_adj < 0.001, "***",
                     ifelse(p_adj < 0.01, "**",
                            ifelse(p_adj < 0.05, "*", "ns")))
  ) %>%
  dplyr::select(comparison, diff, lwr, upr, p_adj, p_signif)

# Display Tukey test results
tank_normalized_tukey_summary <- tank_normalized_tukey_df %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Tukey HSD Test Results - Tank-Level Normalized Worm Burden",
    subtitle = "Post-hoc pairwise comparisons for worms per infected fish"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(diff, lwr, upr, p_adj),
    decimals = 3
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = p_signif,
      rows = p_signif != "ns"
    )
  )
Table

Infection Prevalence

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Summary
Infection Prevalence Contingency Table
Number of infected vs uninfected fish by treatment
Treatment n_infected n_uninfected n_fish percent_infected percent_uninfected
A- T- P+ 31 51 82 37.8 62.2
A+ T- P+ 19 61 80 23.8 76.2
A- T+ P+ 21 56 77 27.3 72.7
A+ T+ P+ 15 55 70 21.4 78.6
Infection Chi-Square
Chi-Square Test Results - Infection Prevalence
Test for independence between treatment and infection status
Statistic Value
Chi-square statistic 6.16510
Degrees of freedom 3.00000
P-value 0.10385
Infection Pairwise
Pairwise Chi-Square Tests - Infection Prevalence
Post-hoc comparisons between treatment pairs
Comparison Chi_square P_value Significant
A- T- P+ vs A+ T- P+ 3.1190 0.077384 No
A- T- P+ vs A- T+ P+ 1.5515 0.212910 No
A- T- P+ vs A+ T+ P+ 4.0541 0.044064 Yes
A+ T- P+ vs A- T+ P+ 0.1045 0.746535 No
A+ T- P+ vs A+ T+ P+ 0.0205 0.886027 No
A- T+ P+ vs A+ T+ P+ 0.3980 0.528098 No

Worm Burden

Summary
Worm Burden Summary by Treatment
Total, mean, median, standard deviation, and standard error of worm counts per treatment
Treatment total_worms mean_worms median_worms sd_worms se_worms n_fish
A- T- P+ 167 2.04 0.00 5.16 0.57 82
A+ T- P+ 45 0.56 0.00 1.58 0.18 80
A- T+ P+ 121 1.57 0.00 3.59 0.41 77
A+ T+ P+ 77 1.10 0.00 2.82 0.34 70
Negative Binomial Model
Negative Binomial Regression Results - Worm Burden
Model: glm.nb(worm counts ~ Treatment)
Treatment Estimate Std_Error z_value p_value p_signif
Intercept 0.711 0.292 2.433 0.015 *
A+ T- P+ −1.287 0.435 −2.958 0.003 **
A- T+ P+ −0.259 0.422 −0.614 0.539 ns
A+ T+ P+ −0.616 0.438 −1.407 0.159 ns
Model Fit Stats
Negative Binomial Model Fit Statistics
Model diagnostics and fit measures
Statistic Value
Log-likelihood -388.0400
AIC 786.0800
Theta 0.1534
SE Theta 0.0226
Deviance 189.0100
Residual Deviance 189.0100
Worm Burden Pairwise
Pairwise Comparisons - Worm Burden (Negative Binomial)
Bonferroni-adjusted pairwise comparisons between treatments
contrast estimate SE z.ratio p.value p_signif
(A- T- P+) - (A+ T- P+) 1.287 0.435 2.958 0.019 *
(A- T- P+) - (A- T+ P+) 0.259 0.422 0.614 1.000 ns
(A- T- P+) - (A+ T+ P+) 0.616 0.438 1.407 0.956 ns
(A+ T- P+) - (A- T+ P+) −1.027 0.443 −2.317 0.123 ns
(A+ T- P+) - (A+ T+ P+) −0.671 0.458 −1.464 0.859 ns
(A- T+ P+) - (A+ T+ P+) 0.357 0.446 0.800 1.000 ns
Model Diagnostics
Negative Binomial Model Diagnostics
Checking model adequacy and dispersion
Statistic Value
Theta (dispersion parameter) 0.1534
SE Theta 0.0226
Deviance 189.01
Degrees of Freedom 305
Deviance/DF 0.62
Model adequacy Good fit
Model Comparison
Model Comparison: Poisson vs Negative Binomial
Comparing model fit and adequacy
Model Log_Likelihood AIC Deviance DF Deviance_DF_Ratio
Poisson −800.23 1,608.45 1,338.72 305 4.39
Negative Binomial −388.04 786.08 189.01 305 0.62

Normalized Worm Burden

Summary
Normalized Worm Burden Summary by Treatment
Worms per infected fish (normalized for infection prevalence)
Treatment total_worms n_infected worms_per_infected n_fish infection_rate
A- T- P+ 167 31 5.39 82 37.80
A+ T- P+ 45 19 2.37 80 23.80
A- T+ P+ 121 21 5.76 77 27.30
A+ T+ P+ 77 15 5.13 70 21.40
ANOVA Results
ANOVA Results - Normalized Worm Burden
One-way ANOVA: worms per infected fish ~ Treatment (tank-level data)
Statistic Value
F-value 1.857900
Degrees of freedom (Treatment) 3.000000
Degrees of freedom (Residuals) 8.000000
P-value 0.215039
Tukey HSD Test
Tukey HSD Test Results - Normalized Worm Burden
Post-hoc pairwise comparisons for worms per infected fish (tank-level data)
comparison diff lwr upr p_adj p_signif
A+ T- P+-A- T- P+ −3.012 −8.471 2.448 0.354 ns
A- T+ P+-A- T- P+ 0.198 −5.262 5.657 0.999 ns
A+ T+ P+-A- T- P+ 0.536 −4.924 5.995 0.988 ns
A- T+ P+-A+ T- P+ 3.210 −2.250 8.669 0.307 ns
A+ T+ P+-A+ T- P+ 3.548 −1.912 9.007 0.238 ns
A+ T+ P+-A- T+ P+ 0.338 −5.122 5.798 0.997 ns
Kruskal-Wallis Test
Kruskal-Wallis Test Results - Normalized Worm Burden
Non-parametric test for differences in worms per infected fish across treatments (tank-level data)
Statistic Value
Kruskal-Wallis chi-squared 3.769200
Degrees of freedom 3.000000
P-value 0.287485
Permutation Test
Permutation Test Results - Normalized Worm Burden
Non-parametric permutation test for ratio data (tank-level data)
Statistic Value
Z statistic 4.5168
P-value (approximate) 0.2356
Number of permutations 10000.0000
Bootstrap Confidence Intervals
Bootstrap Confidence Intervals - Normalized Worm Burden
95% confidence intervals for worms per infected fish by treatment (tank-level data)
Treatment Mean_Ratio Bootstrap_Mean CI_Lower CI_Upper N_Tanks
A- T- P+ 5.095 5.095 3.286 6.167 3
A+ T- P+ 2.083 2.083 1.000 3.625 3
A- T+ P+ 5.293 5.293 2.333 7.545 3
A+ T+ P+ 5.631 5.631 3.143 8.000 3
Pairwise Permutation Tests
Pairwise Permutation Tests - Normalized Worm Burden
Non-parametric pairwise comparisons for ratio data (tank-level data)
Comparison Z_statistic P_value Significant
A- T- P+ vs A+ T- P+ 1.7453 0.2013 No
A- T- P+ vs A- T+ P+ -0.1230 0.9007 No
A- T- P+ vs A+ T+ P+ -0.3536 1.0000 No
A+ T- P+ vs A- T+ P+ -1.5176 0.2005 No
A+ T- P+ vs A+ T+ P+ -1.6551 0.2003 No
A- T+ P+ vs A+ T+ P+ -0.1804 0.8006 No
Model Comparison
Model Comparison - Normalized Worm Burden Analysis
Comparing different statistical approaches for ratio data
Test Statistic P_value Assumptions
ANOVA F = 1.8579 0.215039 Normality, homogeneity of variance
Kruskal-Wallis χ² = 3.7692 0.287485 None (non-parametric)
Permutation Test Z = 4.5168 0.235600 None (non-parametric)
Tank-Level Analysis
Tank-Level Normalized Worm Burden Summary
Worms per infected fish by treatment and tank
Treatment Tank.ID total_worms n_infected worms_per_infected n_fish infection_rate
A- T- P+ 4 70 12 5.83 29 41.40
A- T- P+ 5 23 7 3.29 25 28.00
A- T- P+ 6 74 12 6.17 28 42.90
A+ T- P+ 16 13 8 1.62 29 27.60
A+ T- P+ 17 29 8 3.62 28 28.60
A+ T- P+ 18 3 3 1.00 23 13.00
A- T+ P+ 10 24 4 6.00 21 19.00
A- T+ P+ 11 14 6 2.33 30 20.00
A- T+ P+ 12 83 11 7.55 26 42.30
A+ T+ P+ 22 22 7 3.14 28 25.00
A+ T+ P+ 23 23 4 5.75 23 17.40
A+ T+ P+ 24 32 4 8.00 19 21.10
ANOVA Results - Tank-Level Normalized Worm Burden
One-way ANOVA: worms per infected fish ~ Treatment
Statistic Value
F-value 1.857900
Degrees of freedom (Treatment) 3.000000
Degrees of freedom (Residuals) 8.000000
P-value 0.215039
Tukey HSD Test Results - Tank-Level Normalized Worm Burden
Post-hoc pairwise comparisons for worms per infected fish
comparison diff lwr upr p_adj p_signif
A+ T- P+-A- T- P+ −3.012 −8.471 2.448 0.354 ns
A- T+ P+-A- T- P+ 0.198 −5.262 5.657 0.999 ns
A+ T+ P+-A- T- P+ 0.536 −4.924 5.995 0.988 ns
A- T+ P+-A+ T- P+ 3.210 −2.250 8.669 0.307 ns
A+ T+ P+-A+ T- P+ 3.548 −1.912 9.007 0.238 ns
A+ T+ P+-A- T+ P+ 0.338 −5.122 5.798 0.997 ns

03) Is host-microbiome recovery historically contingent?

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality

All

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

By Tank

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

Statistical Analysis - Recovery Mortality

Code
# Prepare data for statistical analysis - recovery mortality
recovery_mortality_data <- data.mortality %>%
  dplyr::filter(Treatment %in% c("A- T- P-", "A+ T- P-", "A- T+ P-", "A+ T+ P-")) %>%
  dplyr::group_by(Treatment) %>%
  dplyr::count() %>%
  dplyr::mutate(
    total_per_group = 90,
    dead = total_per_group - n,
    percent_mortality = round((dead / total_per_group) * 100, 1)
  )

# Chi-square test for independence between treatment and mortality
set.seed(42)
recovery_chi_square_test <- chisq.test(
  matrix(c(recovery_mortality_data$dead, recovery_mortality_data$n), 
         ncol = 2, byrow = FALSE)
)

# Create mortality rate summary table for display
recovery_mortality_rate_table <- recovery_mortality_data %>%
  dplyr::select(Treatment, dead, n, percent_mortality) %>%
  dplyr::rename(Survived = n) %>%
  dplyr::mutate(
    Total = dead + Survived,
    percent_survived = round((Survived / Total) * 100, 1)
  ) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Mortality Rate Summary - Recovery",
    subtitle = "Percent mortality and survival rates by treatment (parasite-unexposed only)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(dead, Survived, Total),
    decimals = 0
  ) %>%
  gt::fmt_number(
    columns = c(percent_mortality, percent_survived),
    decimals = 1
  )

# Chi-square test results
recovery_chi_square_results <- data.frame(
  Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
  Value = c(
    round(recovery_chi_square_test$statistic, 4),
    recovery_chi_square_test$parameter,
    round(recovery_chi_square_test$p.value, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Chi-Square Test Results - Recovery Mortality",
    subtitle = "Test for independence between treatment and mortality (parasite-unexposed only)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Post-hoc pairwise comparisons using chi-square tests
set.seed(42)
recovery_pairwise_results <- list()
recovery_treatments <- recovery_mortality_data$Treatment

for(i in 1:(length(recovery_treatments)-1)) {
  for(j in (i+1):length(recovery_treatments)) {
    trt1 <- recovery_treatments[i]
    trt2 <- recovery_treatments[j]
    
    # Get data for these two treatments
    data1 <- recovery_mortality_data %>% dplyr::filter(Treatment == trt1)
    data2 <- recovery_mortality_data %>% dplyr::filter(Treatment == trt2)
    
    # Create 2x2 contingency table
    cont_table <- matrix(
      c(data1$dead, data1$n, data2$dead, data2$n),
      nrow = 2, ncol = 2,
      dimnames = list(
        c("Dead", "Survived"),
        c(trt1, trt2)
      )
    )
    
    # Perform chi-square test
    test_result <- chisq.test(cont_table)
    
    recovery_pairwise_results[[paste(trt1, "vs", trt2)]] <- data.frame(
      Comparison = paste(trt1, "vs", trt2),
      Chi_square = round(test_result$statistic, 4),
      P_value = round(test_result$p.value, 6),
      Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
    )
  }
}

# Combine pairwise results
recovery_pairwise_summary <- do.call(rbind, recovery_pairwise_results) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pairwise Chi-Square Tests - Recovery Mortality",
    subtitle = "Post-hoc comparisons between treatment pairs (parasite-unexposed only)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )
Table

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality Rates
Mortality Rate Summary - Recovery
Percent mortality and survival rates by treatment (parasite-unexposed only)
dead Survived percent_mortality Total percent_survived
A- T- P-
9 81 10.0 90 90.0
A+ T- P-
14 76 15.6 90 84.4
A- T+ P-
19 71 21.1 90 78.9
A+ T+ P-
30 60 33.3 90 66.7
Chi-Square
Chi-Square Test Results - Recovery Mortality
Test for independence between treatment and mortality (parasite-unexposed only)
Statistic Value
Chi-square statistic 16.805600
Degrees of freedom 3.000000
P-value 0.000775
Chi-Square Pairwise
Pairwise Chi-Square Tests - Recovery Mortality
Post-hoc comparisons between treatment pairs (parasite-unexposed only)
Comparison Chi_square P_value Significant
A- T- P- vs A+ T- P- 0.7976 0.371823 No
A- T- P- vs A- T+ P- 3.4258 0.064187 No
A- T- P- vs A+ T+ P- 13.0933 0.000296 Yes
A+ T- P- vs A- T+ P- 0.5937 0.440995 No
A+ T- P- vs A+ T+ P- 6.7680 0.009280 Yes
A- T+ P- vs A+ T+ P- 2.8042 0.094019 No

04) History of stressors

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality

All

By Tank

Statistical Analysis - History of Stressors Mortality

Code
# Prepare data for statistical analysis - percent mortality rates by History
history_mortality_stats_data <- data.mortality %>%
  dplyr::group_by(History) %>%
  dplyr::count() %>%
  dplyr::mutate(
    # Calculate expected totals: 90 fish per treatment group
    total_per_group = dplyr::case_when(
      History == 0 ~ 90 * 2,  # 2 treatment groups: A- T- P- and A- T- P+
      History == 1 ~ 90 * 4,  # 4 treatment groups: A+ T- P-, A+ T- P+, A- T+ P-, A- T+ P+
      History == 2 ~ 90 * 2,  # 2 treatment groups: A+ T+ P- and A+ T+ P+
      TRUE ~ 0
    ),
    dead = total_per_group - n,
    percent_mortality = round((dead / total_per_group) * 100, 1),
    History_Label = dplyr::case_when(
      History == 0 ~ "No prior stressors",
      History == 1 ~ "One prior stressor", 
      History == 2 ~ "Two prior stressors",
      TRUE ~ "Unknown"
    )
  )

# Chi-square test for independence between history and mortality
set.seed(42)
history_chi_square_test <- chisq.test(
  matrix(c(history_mortality_stats_data$dead, history_mortality_stats_data$n), 
         ncol = 2, byrow = FALSE)
)

# Create mortality rate summary table for display
history_mortality_rate_table <- history_mortality_stats_data %>%
  dplyr::select(History, History_Label, dead, n, percent_mortality) %>%
  dplyr::rename(Survived = n) %>%
  dplyr::mutate(
    Total = dead + Survived,
    percent_survived = round((Survived / Total) * 100, 1)
  ) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Mortality Rate Summary - History of Stressors",
    subtitle = "Percent mortality and survival rates by history of stressors"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(dead, Survived, Total),
    decimals = 0
  ) %>%
  gt::fmt_number(
    columns = c(percent_mortality, percent_survived),
    decimals = 1
  )

# Chi-square test results
history_test_stat_table <- data.frame(
  Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
  Value = c(
    round(history_chi_square_test$statistic, 4),
    history_chi_square_test$parameter,
    round(history_chi_square_test$p.value, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Chi-Square Test Results - History of Stressors Mortality",
    subtitle = "Test for independence between history of stressors and mortality"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Post-hoc pairwise comparisons using chi-square tests
set.seed(42)
history_pairwise_results <- list()
history_levels <- history_mortality_stats_data$History

for(i in 1:(length(history_levels)-1)) {
  for(j in (i+1):length(history_levels)) {
    hist1 <- history_levels[i]
    hist2 <- history_levels[j]
    
    # Get data for these two history levels
    data1 <- history_mortality_stats_data %>% dplyr::filter(History == hist1)
    data2 <- history_mortality_stats_data %>% dplyr::filter(History == hist2)
    
    # Create 2x2 contingency table
    cont_table <- matrix(
      c(data1$dead, data1$n, data2$dead, data2$n),
      nrow = 2, ncol = 2,
      dimnames = list(
        c("Dead", "Survived"),
        c(paste("History", hist1), paste("History", hist2))
      )
    )
    
    # Perform chi-square test
    test_result <- chisq.test(cont_table)
    
    history_pairwise_results[[paste("History", hist1, "vs", "History", hist2)]] <- data.frame(
      Comparison = paste("History", hist1, "vs", "History", hist2),
      Chi_square = round(test_result$statistic, 4),
      P_value = round(test_result$p.value, 6),
      Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
    )
  }
}

# Combine pairwise results
history_pairwise_summary <- do.call(rbind, history_pairwise_results) %>%
  dplyr::arrange(desc(Significant), P_value) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pairwise Chi-Square Tests - History of Stressors Mortality",
    subtitle = "Post-hoc comparisons between history levels"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )

# Show only significant pairwise comparisons
history_pairwise_summary__SigOnly <- do.call(rbind, history_pairwise_results) %>%
  dplyr::arrange(desc(Significant), P_value) %>%
  dplyr::filter(Significant == "Yes") %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Significant Pairwise Chi-Square Tests - History of Stressors Mortality",
    subtitle = "Post-hoc comparisons between history levels (significant only)"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )

# Additional analysis: Trend analysis across history levels
# Calculate trend in mortality rates
history_trend_data <- history_mortality_stats_data %>%
  dplyr::arrange(History) %>%
  dplyr::select(History, percent_mortality)

# Linear trend test (Cochran-Armitage trend test equivalent using chi-square)
set.seed(42)
# Create trend weights (0, 1, 2 for history levels)
trend_weights <- c(0, 1, 2)
trend_chi_square <- chisq.test(
  matrix(c(history_mortality_stats_data$dead, history_mortality_stats_data$n), 
         ncol = 2, byrow = FALSE),
  simulate.p.value = TRUE
)

# Calculate correlation between history level and mortality rate
history_correlation <- cor.test(
  history_trend_data$History, 
  history_trend_data$percent_mortality, 
  method = "pearson"
)

# Trend analysis results
history_trend_results <- data.frame(
  Statistic = c("Correlation coefficient (r)", "P-value (correlation)", "Chi-square trend test p-value"),
  Value = c(
    round(history_correlation$estimate, 4),
    round(history_correlation$p.value, 6),
    round(trend_chi_square$p.value, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Trend Analysis - History of Stressors Mortality",
    subtitle = "Testing for linear trend in mortality rates across history levels"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Effect size analysis: Calculate odds ratios between history levels
history_effect_sizes <- data.frame(
  Comparison = c("History 1 vs History 0", "History 2 vs History 0", "History 2 vs History 1"),
  Odds_Ratio = c(
    round((history_mortality_stats_data$dead[2] * history_mortality_stats_data$n[1]) / 
          (history_mortality_stats_data$dead[1] * history_mortality_stats_data$n[2]), 3),
    round((history_mortality_stats_data$dead[3] * history_mortality_stats_data$n[1]) / 
          (history_mortality_stats_data$dead[1] * history_mortality_stats_data$n[3]), 3),
    round((history_mortality_stats_data$dead[3] * history_mortality_stats_data$n[2]) / 
          (history_mortality_stats_data$dead[2] * history_mortality_stats_data$n[3]), 3)
  ),
  Mortality_Rate_1 = c(
    round(history_mortality_stats_data$percent_mortality[2], 1),
    round(history_mortality_stats_data$percent_mortality[3], 1),
    round(history_mortality_stats_data$percent_mortality[3], 1)
  ),
  Mortality_Rate_2 = c(
    round(history_mortality_stats_data$percent_mortality[1], 1),
    round(history_mortality_stats_data$percent_mortality[1], 1),
    round(history_mortality_stats_data$percent_mortality[2], 1)
  )
) %>%
  dplyr::mutate(
    Interpretation = dplyr::case_when(
      Odds_Ratio > 1 ~ "Higher mortality in first group",
      Odds_Ratio < 1 ~ "Lower mortality in first group",
      Odds_Ratio == 1 ~ "No difference"
    )
  ) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Effect Size Analysis - History of Stressors Mortality",
    subtitle = "Odds ratios comparing mortality rates between history levels"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(Odds_Ratio, Mortality_Rate_1, Mortality_Rate_2),
    decimals = 1
  )
Table

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality Rates
Mortality Rate Summary - History of Stressors
Percent mortality and survival rates by history of stressors
History_Label dead Survived percent_mortality Total percent_survived
0
No prior stressors 17 163 9.4 180 90.6
1
One prior stressor 56 304 15.6 360 84.4
2
Two prior stressors 50 130 27.8 180 72.2
Chi-Square
Chi-Square Test Results - History of Stressors Mortality
Test for independence between history of stressors and mortality
Statistic Value
Chi-square statistic 22.542000
Degrees of freedom 2.000000
P-value 0.000013
Chi-Square Pairwise
Pairwise Chi-Square Tests - History of Stressors Mortality
Post-hoc comparisons between history levels
Comparison Chi_square P_value Significant
History 0 vs History 2 18.7785 0.000015 Yes
History 1 vs History 2 10.6010 0.001130 Yes
History 0 vs History 1 3.3284 0.068094 No
Significant Pairwise Chi-Square Tests - History of Stressors Mortality
Post-hoc comparisons between history levels (significant only)
Comparison Chi_square P_value Significant
History 0 vs History 2 18.7785 0.000015 Yes
History 1 vs History 2 10.6010 0.001130 Yes
Trend Analysis
Trend Analysis - History of Stressors Mortality
Testing for linear trend in mortality rates across history levels
Statistic Value
Correlation coefficient (r) 0.982700
P-value (correlation) 0.118467
Chi-square trend test p-value 0.000500
Effect Sizes
Effect Size Analysis - History of Stressors Mortality
Odds ratios comparing mortality rates between history levels
Comparison Odds_Ratio Mortality_Rate_1 Mortality_Rate_2 Interpretation
History 1 vs History 0 1.8 15.6 9.4 Higher mortality in first group
History 2 vs History 0 3.7 27.8 9.4 Higher mortality in first group
History 2 vs History 1 2.1 27.8 15.6 Higher mortality in first group

05) History of stressors by pathogen exposure

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Mortality

All

By Tank

Statistical Analysis - History × Pathogen Interaction

Code
# Prepare data for statistical analysis - History × Pathogen interaction
interaction_stats_data <- history_pathogen_interaction_data %>%
  dplyr::select(History, Pathogen, History_Label, Pathogen_Label, dead, n, percent_mortality) %>%
  dplyr::rename(Survived = n) %>%
  dplyr::mutate(
    Total = dead + Survived,
    percent_survived = round((Survived / Total) * 100, 1)
  )

# 1. Chi-square test for independence between History and Pathogen exposure
set.seed(42)
# Create 3x2 contingency table (History × Pathogen)
interaction_contingency <- matrix(
  c(interaction_stats_data$dead), 
  nrow = 3, ncol = 2,
  dimnames = list(
    c("History 0", "History 1", "History 2"),
    c("Unexposed", "Exposed")
  )
)

interaction_chi_square <- chisq.test(interaction_contingency)

# 2. Test for interaction effect using logistic regression approach
# Create individual-level data for interaction analysis
set.seed(42)
individual_interaction_data <- data.mortality %>%
  dplyr::select(History, Pathogen, Sample) %>%
  dplyr::mutate(
    # Create binary outcome: 1 = survived, 0 = died (assuming all sampled fish survived)
    Outcome = 1,  # All fish in the dataset survived (were sampled)
    History = factor(History, levels = c(0, 1, 2)),
    Pathogen = factor(Pathogen, levels = c(0, 1))
  )

# Fit logistic regression model with interaction
interaction_model <- glm(Outcome ~ History * Pathogen, 
                        data = individual_interaction_data, 
                        family = "binomial")
## Warning: glm.fit: algorithm did not converge
# Get model summary
interaction_summary <- summary(interaction_model)

# Extract coefficients and significance
interaction_coefficients <- data.frame(
  Term = rownames(interaction_summary$coefficients),
  Estimate = round(interaction_summary$coefficients[, 1], 4),
  Std_Error = round(interaction_summary$coefficients[, 2], 4),
  z_value = round(interaction_summary$coefficients[, 3], 4),
  p_value = round(interaction_summary$coefficients[, 4], 6)
) %>%
  dplyr::mutate(
    p_signif = ifelse(p_value < 0.001, "***", 
                     ifelse(p_value < 0.01, "**",
                            ifelse(p_value < 0.05, "*", "ns")))
  )

# 3. Create summary tables for display
interaction_summary_table <- interaction_stats_data %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Mortality Summary - History × Pathogen Interaction",
    subtitle = "Percent mortality and survival rates by history and pathogen exposure"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(dead, Survived, Total),
    decimals = 0
  ) %>%
  gt::fmt_number(
    columns = c(percent_mortality, percent_survived),
    decimals = 1
  )

# Chi-square test results
interaction_chi_results <- data.frame(
  Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
  Value = c(
    round(interaction_chi_square$statistic, 4),
    interaction_chi_square$parameter,
    round(interaction_chi_square$p.value, 6)
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Chi-Square Test Results - History × Pathogen Interaction",
    subtitle = "Test for independence between history and pathogen exposure"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# Logistic regression results
interaction_model_results <- interaction_coefficients %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Logistic Regression Results - History × Pathogen Interaction",
    subtitle = "Model: Outcome ~ History * Pathogen"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(Estimate, Std_Error, z_value, p_value),
    decimals = 3
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = p_signif,
      rows = p_signif != "ns"
    )
  )

# 4. Calculate effect sizes for each comparison
# Calculate odds ratios for each history level comparing exposed vs unexposed
set.seed(42)
effect_size_comparisons <- data.frame(
  Comparison = c(
    "History 0: Exposed vs Unexposed",
    "History 1: Exposed vs Unexposed", 
    "History 2: Exposed vs Unexposed"
  ),
  Odds_Ratio = c(
    # History 0: Exposed vs Unexposed
    round((interaction_stats_data$dead[2] * interaction_stats_data$Survived[1]) / 
          (interaction_stats_data$dead[1] * interaction_stats_data$Survived[2]), 3),
    # History 1: Exposed vs Unexposed
    round((interaction_stats_data$dead[4] * interaction_stats_data$Survived[3]) / 
          (interaction_stats_data$dead[3] * interaction_stats_data$Survived[4]), 3),
    # History 2: Exposed vs Unexposed
    round((interaction_stats_data$dead[6] * interaction_stats_data$Survived[5]) / 
          (interaction_stats_data$dead[5] * interaction_stats_data$Survived[6]), 3)
  ),
  Mortality_Exposed = c(
    round(interaction_stats_data$percent_mortality[2], 1),
    round(interaction_stats_data$percent_mortality[4], 1),
    round(interaction_stats_data$percent_mortality[6], 1)
  ),
  Mortality_Unexposed = c(
    round(interaction_stats_data$percent_mortality[1], 1),
    round(interaction_stats_data$percent_mortality[3], 1),
    round(interaction_stats_data$percent_mortality[5], 1)
  )
) %>%
  dplyr::mutate(
    Mortality_Difference = Mortality_Exposed - Mortality_Unexposed,
    Interpretation = dplyr::case_when(
      Odds_Ratio > 1 ~ "Higher mortality when exposed",
      Odds_Ratio < 1 ~ "Lower mortality when exposed",
      Odds_Ratio == 1 ~ "No difference"
    )
  ) 

effect_size_comparisons_table <- effect_size_comparisons %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Effect Size Analysis - History × Pathogen Interaction",
    subtitle = "Odds ratios comparing exposed vs unexposed within each history level"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(Odds_Ratio, Mortality_Exposed, Mortality_Unexposed, Mortality_Difference),
    decimals = 1
  )

# 5. Test for trend in pathogen effect across history levels
set.seed(42)
# Calculate correlation between history level and the difference in mortality (exposed - unexposed)
trend_data <- effect_size_comparisons %>%
  dplyr::mutate(History_Level = c(0, 1, 2))

pathogen_effect_trend <- cor.test(
  trend_data$History_Level, 
  trend_data$Mortality_Difference, 
  method = "pearson"
)

# Trend analysis results
pathogen_effect_trend_results <- data.frame(
  Statistic = c("Correlation coefficient (r)", "P-value", "Interpretation"),
  Value = c(
    round(pathogen_effect_trend$estimate, 4),
    round(pathogen_effect_trend$p.value, 6),
    ifelse(pathogen_effect_trend$p.value < 0.05, 
           "Significant trend in pathogen effect across history levels",
           "No significant trend in pathogen effect across history levels")
  )
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Trend Analysis - Pathogen Effect Across History Levels",
    subtitle = "Testing whether the effect of pathogen exposure changes across history levels"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  )

# 6. Post-hoc pairwise comparisons within each pathogen group
set.seed(42)
# Within unexposed group
unexposed_data <- interaction_stats_data %>% dplyr::filter(Pathogen == 0)
unexposed_pairwise <- list()

for(i in 1:(nrow(unexposed_data)-1)) {
  for(j in (i+1):nrow(unexposed_data)) {
    hist1 <- unexposed_data$History[i]
    hist2 <- unexposed_data$History[j]
    
    # Create 2x2 contingency table
    cont_table <- matrix(
      c(unexposed_data$dead[i], unexposed_data$Survived[i], 
        unexposed_data$dead[j], unexposed_data$Survived[j]),
      nrow = 2, ncol = 2,
      dimnames = list(
        c("Dead", "Survived"),
        c(paste("History", hist1), paste("History", hist2))
      )
    )
    
    # Perform chi-square test
    test_result <- chisq.test(cont_table)
    
    unexposed_pairwise[[paste("History", hist1, "vs", "History", hist2)]] <- data.frame(
      Comparison = paste("History", hist1, "vs", "History", hist2, "(Unexposed)"),
      Chi_square = round(test_result$statistic, 4),
      P_value = round(test_result$p.value, 6),
      Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
    )
  }
}

# Within exposed group
exposed_data <- interaction_stats_data %>% dplyr::filter(Pathogen == 1)
exposed_pairwise <- list()

for(i in 1:(nrow(exposed_data)-1)) {
  for(j in (i+1):nrow(exposed_data)) {
    hist1 <- exposed_data$History[i]
    hist2 <- exposed_data$History[j]
    
    # Create 2x2 contingency table
    cont_table <- matrix(
      c(exposed_data$dead[i], exposed_data$Survived[i], 
        exposed_data$dead[j], exposed_data$Survived[j]),
      nrow = 2, ncol = 2,
      dimnames = list(
        c("Dead", "Survived"),
        c(paste("History", hist1), paste("History", hist2))
      )
    )
    
    # Perform chi-square test
    test_result <- chisq.test(cont_table)
    
    exposed_pairwise[[paste("History", hist1, "vs", "History", hist2)]] <- data.frame(
      Comparison = paste("History", hist1, "vs", "History", hist2, "(Exposed)"),
      Chi_square = round(test_result$statistic, 4),
      P_value = round(test_result$p.value, 6),
      Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
    )
  }
}

# 7. Pairwise comparisons between same history types across pathogen exposures
set.seed(42)
same_history_pathogen_pairwise <- list()

# For each history level, compare unexposed vs exposed
for(hist_level in c(0, 1, 2)) {
  # Get data for this history level
  unexposed_row <- interaction_stats_data %>% dplyr::filter(History == hist_level & Pathogen == 0)
  exposed_row <- interaction_stats_data %>% dplyr::filter(History == hist_level & Pathogen == 1)
  
  if(nrow(unexposed_row) > 0 && nrow(exposed_row) > 0) {
    # Calculate percent mortality for each group
    unexposed_mortality <- unexposed_row$percent_mortality
    exposed_mortality <- exposed_row$percent_mortality
    
    # Calculate the difference in percent mortality
    mortality_difference <- exposed_mortality - unexposed_mortality
    
    # Perform t-test on percent mortality rates
    # Since we have single values, we'll use a simple comparison
    # For small sample sizes, we can use a z-test for proportions
    # Calculate standard error of the difference
    unexposed_se <- sqrt((unexposed_mortality * (100 - unexposed_mortality)) / unexposed_row$Total)
    exposed_se <- sqrt((exposed_mortality * (100 - exposed_mortality)) / exposed_row$Total)
    se_diff <- sqrt(unexposed_se^2 + exposed_se^2)
    
    # Calculate z-statistic
    z_stat <- mortality_difference / se_diff
    p_value <- 2 * (1 - pnorm(abs(z_stat)))
    
    same_history_pathogen_pairwise[[paste("History", hist_level)]] <- data.frame(
      Comparison = paste("History", hist_level, "(Unexposed) vs History", hist_level, "(Exposed)"),
      Unexposed_Mortality = round(unexposed_mortality, 1),
      Exposed_Mortality = round(exposed_mortality, 1),
      Mortality_Difference = round(mortality_difference, 1),
      Z_statistic = round(z_stat, 4),
      P_value = round(p_value, 6),
      Significant = ifelse(p_value < 0.05, "Yes", "No")
    )
  }
}

# Create separate tables for different types of pairwise comparisons

# Table 1: Within unexposed group comparisons
unexposed_pairwise_results <- do.call(rbind, unexposed_pairwise) %>%
  dplyr::arrange(desc(Significant), P_value) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pairwise Chi-Square Tests - Within Unexposed Group",
    subtitle = "Comparing history levels within unexposed fish"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )

# Table 2: Within exposed group comparisons
exposed_pairwise_results <- do.call(rbind, exposed_pairwise) %>%
  dplyr::arrange(desc(Significant), P_value) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pairwise Chi-Square Tests - Within Exposed Group",
    subtitle = "Comparing history levels within exposed fish"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )

# Create separate table for same history comparisons across pathogen exposures
same_history_pathogen_results <- do.call(rbind, same_history_pathogen_pairwise) %>%
  dplyr::arrange(desc(Significant), P_value) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Pathogen Effect Within Each History Level",
    subtitle = "Comparing percent mortality between unexposed vs exposed within the same history of stressors"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "#f7f7f7"),
    locations = gt::cells_column_labels()
  ) %>%
  gt::fmt_number(
    columns = c(Unexposed_Mortality, Exposed_Mortality, Mortality_Difference),
    decimals = 1
  ) %>%
  gt::fmt_number(
    columns = c(Z_statistic, P_value),
    decimals = 4
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = Significant,
      rows = Significant == "Yes"
    )
  )
Table

(hide)

Click on tabs to display tables. Scroll to see additional rows.

Summary Table
Mortality Summary - History × Pathogen Interaction
Percent mortality and survival rates by history and pathogen exposure
History_Label Pathogen_Label dead Survived percent_mortality Total percent_survived
0 - 0
No prior stressors Unexposed 9 81 10.0 90 90.0
0 - 1
No prior stressors Exposed 8 82 8.9 90 91.1
1 - 0
One prior stressor Unexposed 33 147 18.3 180 81.7
1 - 1
One prior stressor Exposed 23 157 12.8 180 87.2
2 - 0
Two prior stressors Unexposed 30 60 33.3 90 66.7
2 - 1
Two prior stressors Exposed 20 70 22.2 90 77.8
Chi-Square Test
Chi-Square Test Results - History × Pathogen Interaction
Test for independence between history and pathogen exposure
Statistic Value
Chi-square statistic 18.392800
Degrees of freedom 2.000000
P-value 0.000101
Logistic Regression
Logistic Regression Results - History × Pathogen Interaction
Model: Outcome ~ History * Pathogen
Term Estimate Std_Error z_value p_value p_signif
(Intercept) 26.566 39,569.344 0.001 0.999 ns
History1 0.000 49,279.608 0.000 1.000 ns
History2 0.000 60,658.574 0.000 1.000 ns
Pathogen1 0.000 55,788.569 0.000 1.000 ns
History1:Pathogen1 0.000 69,158.557 0.000 1.000 ns
History2:Pathogen1 0.000 83,891.969 0.000 1.000 ns
Effect Sizes
Effect Size Analysis - History × Pathogen Interaction
Odds ratios comparing exposed vs unexposed within each history level
Comparison Odds_Ratio Mortality_Exposed Mortality_Unexposed Mortality_Difference Interpretation
History 0: Exposed vs Unexposed 0.9 8.9 10.0 −1.1 Lower mortality when exposed
History 1: Exposed vs Unexposed 0.7 12.8 18.3 −5.5 Lower mortality when exposed
History 2: Exposed vs Unexposed 0.6 22.2 33.3 −11.1 Lower mortality when exposed
Trend Analysis
Trend Analysis - Pathogen Effect Across History Levels
Testing whether the effect of pathogen exposure changes across history levels
Statistic Value
Correlation coefficient (r) -0.9976
P-value 0.044036
Interpretation Significant trend in pathogen effect across history levels
Within Unexposed Group
Pairwise Chi-Square Tests - Within Unexposed Group
Comparing history levels within unexposed fish
Comparison Chi_square P_value Significant
History 0 vs History 2 (Unexposed) 13.0933 0.000296 Yes
History 1 vs History 2 (Unexposed) 6.7314 0.009473 Yes
History 0 vs History 1 (Unexposed) 2.5693 0.108955 No
Within Exposed Group
Pairwise Chi-Square Tests - Within Exposed Group
Comparing history levels within exposed fish
Comparison Chi_square P_value Significant
History 0 vs History 2 (Exposed) 5.1175 0.023686 Yes
History 1 vs History 2 (Exposed) 3.3228 0.068326 No
History 0 vs History 1 (Exposed) 0.5512 0.457833 No
Same History Across Pathogen Exposure
Pathogen Effect Within Each History Level
Comparing percent mortality between unexposed vs exposed within the same history of stressors
Comparison Unexposed_Mortality Exposed_Mortality Mortality_Difference Z_statistic P_value Significant
History 2 (Unexposed) vs History 2 (Exposed) 33.3 22.2 −11.1 −1.6759 0.0938 No
History 1 (Unexposed) vs History 1 (Exposed) 18.3 12.8 −5.5 −1.4440 0.1487 No
History 0 (Unexposed) vs History 0 (Exposed) 10.0 8.9 −1.1 −0.2523 0.8008 No